
#include "MCRadiationHeatFluxUserObject.h"
#include "FEProblem.h"
#include "libmesh/elem.h"
#include "libmesh/fe.h"
#include "libmesh/string_to_enum.h"
#include "libmesh/quadrature_gauss.h"
#include "libmesh/plane.h"
#include "MooseMesh.h"
#include <fstream>
#include <time.h>
#include "RayLine.h"
#include "MooseRandom.h"

//#include "MemData.h"
using namespace std;

registerMooseObject("RoshanApp", MCRadiationHeatFluxUserObject);
template<>
InputParameters validParams<MCRadiationHeatFluxUserObject>()
{
	InputParameters params = validParams<SideUserObject>();
	params += validParams<RandomInterface>();
	params.addCoupledVar("temperature", "温度场");
	params.addParam< string >("boundary_groups", "The list of boundary groups");
	params.addParam< bool >("have_medium", false, "have medium or not");
	params.addParam< std::vector<SubdomainName> >("block", "The list of boundary IDs from the mesh where this boundary condition applies");
	params.addParam< int > ("particle_count", 1000000, "单元发射粒子数");
	params.addParam< int >("max_reflect_count", 10, "最大反射次数");
	params.addParam< vector<Real> > ("transmissivity", "透射率");
	params.addParam< vector<Real> > ("absorptivity", "吸收率");
	params.addParam< vector<Real> > ("diffuse_reflectivity", "漫反射百分比");
	params.addParam< vector<Real> > ("mirrors_reflectivity", "镜反射百分比");
	params.addParam< bool >("mpi_or_omp", true, "The style of parallel");
	params.addParam< int >("n_threads", 1, "Number of thread while using openMP");
	params.addParam< int >("BSMC", 0, "Deal with BSMC or not");
	params.addParam< bool >("write_or_read", true,  "Write a file of RD, or read RD from a file");
	params.addParam< string >("filename", "The name of file_RD");
	
	return params;
}

MCRadiationHeatFluxUserObject::MCRadiationHeatFluxUserObject(const InputParameters & parameters) :
	SideUserObject(parameters),
	RandomInterface(parameters, *parameters.get<FEProblem *>("_fe_problem"), parameters.get<THREAD_ID>("_tid"), false),
	_temperature(coupledValue("temperature")),
	_have_medium(getParam< bool >("have_medium")),
	_particle_count(getParam<int> ("particle_count")),
	_max_reflect_count(getParam<int>("max_reflect_count")),
	_transmissivity(getParam< vector<Real> > ("transmissivity")),
	_absorptivity(getParam< vector<Real> > ("absorptivity")),
	_diffuse_reflectivity(getParam< vector<Real> > ("diffuse_reflectivity")),
	_mirrors_reflectivity(getParam< vector<Real> > ("mirrors_reflectivity")),
	_mpi_or_omp(getParam< bool >("mpi_or_omp")),
	_n_threads(getParam< int >("n_threads")),
	_BSMC(getParam< int >("BSMC")),
	_write_or_read(getParam< bool > ("write_or_read")),
	_filename(getParam< string >("filename"))
{
	if(_have_medium)
	{
		vector<SubdomainName> block = getParam<std::vector<SubdomainName> >("block");
		for(vector<SubdomainName>::iterator it = block.begin(); it != block.end(); ++it)
		{
			SubdomainID id = _mesh.getSubdomainID(*it);
			_block_ids.insert(id);
		}
	}
}

MCRadiationHeatFluxUserObject::~MCRadiationHeatFluxUserObject()
{
	for(unsigned int i=0; i<_all_element.size(); i++)
	{
		delete _all_element[i]->_elem;
	}
}

void MCRadiationHeatFluxUserObject::initialSetup()
{
	string boundary_groups = getParam<string>("boundary_groups");
	vector<BoundaryName> group;
	vector< vector<BoundaryName> > group_set;
	MooseUtils::tokenize<BoundaryName>(boundary_groups, group, 1, ",");
	group_set.resize(group.size());
	for(unsigned int i = 0; i < group.size(); ++i)
	{
		MooseUtils::tokenize<BoundaryName>(group[i], group_set[i], 1, " ");
	}

	for(unsigned int i = 0; i < group.size(); ++i)
	{
		cout << "boundary_group_" << i << ": " << group[i] << endl;
	}

	cout << endl << "Mesh modifying..." << endl;

	vector<BoundaryName> boundary = getParam<std::vector<BoundaryName> >("boundary");
	std::set<BoundaryID> boundary_ids;
	for(vector<BoundaryName>::iterator it = boundary.begin(); it != boundary.end(); ++it)
	{
		BoundaryID id = _mesh.getBoundaryID(*it);
		boundary_ids.insert(id);
	}

	MeshBase & mesh = _mesh.getMesh();
	UserDefinedMesh  mymesh(_mesh);
	const BoundaryInfo &bnd_info = mesh.get_boundary_info();

	if(_have_medium)
	{
		int count_sideelement = 0;
		int meshelems = mesh.n_elem();
		
		map<dof_id_type, int> transform;
		int myelemcount = 0;
		for (int nelem = 0; nelem < meshelems; nelem++)
		{
			const Elem * elem = mesh.elem_ptr(nelem);
			SubdomainID curr_subdomain = elem->subdomain_id();
			if( ElemInBlock(curr_subdomain,_block_ids))
			{
				UserDefinedElem newelem;
				newelem._elem = elem;
				transform[elem->id()] = myelemcount;
				myelemcount++;
				mymesh._userDefinedElem.push_back(newelem);
			}
		}
		
		int nelems = mymesh._userDefinedElem.size();
		for (int nelem =0; nelem<nelems; nelem++)
		{
			UserDefinedElem * myelem = &(mymesh._userDefinedElem[nelem]);
			const Elem * elem = myelem->_elem;
			int nsides = elem->n_sides();
			
			for (int nside=0; nside < nsides; nside++)
			{
				UserDefinedSideElem newsideelem;
				myelem->_userDefinedSideElem.push_back(newsideelem);
			}
			
			myelem->_haveWhichSideElement.resize(nsides, -1);
			for (int nside = 0; nside < nsides; nside++)
			{
				UserDefinedSideElem * mysideelem = &(mymesh._userDefinedElem[nelem]._userDefinedSideElem[nside]);
				mysideelem->_elem = elem->side_ptr(nside).release();//.get();
				mysideelem->_left_element = &(mymesh._userDefinedElem[nelem]);

				if( ElemHaveNeighborInBlock(elem->neighbor_ptr(nside), _block_ids) )
				{
					mysideelem->_right_element = &(mymesh._userDefinedElem[ transform[elem->neighbor_ptr(nside)->id()] ]);
				}
				else
				{
					mysideelem->_right_element = NULL;

					const Elem * elem_side = elem->build_side_ptr(nside).release();//.get();
					
					vector<boundary_id_type> bnd_id_vec;
					bnd_info.boundary_ids(elem, nside, bnd_id_vec);
					int bnd_id = bnd_id_vec[0];
					//int bnd_id = bnd_info.boundary_id(elem, nside);

					int bnd_in_which_group =-1;
					for(unsigned int i = 0; i < group.size(); ++i)
					{
						for(unsigned int j = 0; j < group_set[i].size(); ++j)
						{
							if(_mesh.getBoundaryID(group_set[i][j]) == bnd_id)
							{
								bnd_in_which_group = i;
								break;
							}
						}
					}

					if(find(boundary_ids.begin(), boundary_ids.end(), bnd_id) == boundary_ids.end())
					{
						delete elem_side;
						continue;
					}

					unsigned int dim = _mesh.dimension();
					FEType fe_type(Utility::string_to_enum<Order>("CONSTANT"), Utility::string_to_enum<FEFamily>("MONOMIAL"));
					FEBase * _fe_face = (FEBase::build(dim, fe_type)).release();
					QGauss * _qface = new QGauss(dim-1, FIRST);
					_fe_face->attach_quadrature_rule(_qface);
					_fe_face->reinit(elem, nside);
					const std::vector<Point> normals = _fe_face->get_normals();
					_all_element.push_back(new SideElement(elem_side, -normals[0], _transmissivity[bnd_in_which_group], _absorptivity[bnd_in_which_group], _diffuse_reflectivity[bnd_in_which_group], _mirrors_reflectivity[bnd_in_which_group]));
					_all_element[count_sideelement]->_belong_to_which_elem = myelem;
					_all_element[count_sideelement]->_is_which_sideelem = nside;
					myelem->_haveWhichSideElement[nside] = count_sideelement;
					count_sideelement++;

					delete _fe_face;
					delete _qface;

				}
			}
		}
	}
	else
	{
		MeshBase::const_element_iterator el = mesh.active_elements_begin();
		const MeshBase::const_element_iterator end_el = mesh.active_elements_end();
		for ( ; el != end_el ; ++el)
		{
			const Elem *elem = *el;

			for (unsigned int side=0; side < elem->n_sides(); ++side)
			{
				if (elem->neighbor_ptr(side))
					continue;

				const Elem *elem_side = elem->build_side_ptr(side).release();
				vector<boundary_id_type> bnd_id_vec;
				bnd_info.boundary_ids(elem, side, bnd_id_vec);
				int bnd_id = bnd_id_vec[0];
				//int bnd_id = bnd_info.boundary_id(elem, side);

				int bnd_in_which_group =-1;
				for(unsigned int i = 0; i < group.size(); ++i)
				{
					for(unsigned int j = 0; j < group_set[i].size(); ++j)
					{
						if(_mesh.getBoundaryID(group_set[i][j]) == bnd_id)
						{
							bnd_in_which_group = i;
							break;
						}
					}
				}

				if(find(boundary_ids.begin(), boundary_ids.end(), bnd_id) == boundary_ids.end())
				{
					delete elem_side;
					continue;
				}

				unsigned int dim = _mesh.dimension();
				FEType fe_type(Utility::string_to_enum<Order>("CONSTANT"), Utility::string_to_enum<FEFamily>("MONOMIAL"));
				FEBase * _fe_face = (FEBase::build(dim, fe_type)).release();
				QGauss * _qface = new QGauss(dim-1, FIRST);
				_fe_face->attach_quadrature_rule(_qface);
				_fe_face->reinit(elem, side);
				const std::vector<Point> normals = _fe_face->get_normals();
				_all_element.push_back(new SideElement(elem_side, normals[0], _transmissivity[bnd_in_which_group], _absorptivity[bnd_in_which_group], _diffuse_reflectivity[bnd_in_which_group], _mirrors_reflectivity[bnd_in_which_group]));

				delete _fe_face;
				delete _qface;
			}
		}
	}

	for (unsigned int i = 0; i < _all_element.size(); i++)
	{
		_center_point_list.push_back(_all_element[i]->_elem->centroid() );
		_epsilon_list.push_back(_all_element[i]->_absorptivity );
	}

	cout << "Mesh modify completed" << endl;

	cout << endl << "Count of boundary element: " << _all_element.size() << endl << endl;

	_temperature_pow4_bar.resize(_all_element.size(), 0);
	_temperature_pow3_bar.resize(_all_element.size(), 0);
	_temperature_bar.resize(_all_element.size(), 0);
	_flux_radiation.resize(_all_element.size(), 0);
	_flux_radiation_jacobi.resize(_all_element.size(), 0);

	clock_t start_time,end_time;
	start_time = clock();

	if (_write_or_read)
	{
		_communicator.barrier();

		for (unsigned int ii=0;ii<_all_element.size();ii++)
		{
			_all_element[ii]->_local_RD.resize(_all_element.size(), 0);
		}

		if(_mpi_or_omp)
		{
			int local_p;
			int quotient = _particle_count/n_processors();
			int remainder = _particle_count%n_processors();
			if(processor_id() < remainder)
			{
				local_p = quotient + 1;
			}
			else
			{
				local_p = quotient;
			}

			computeRD(local_p, _have_medium);

			_communicator.barrier();

			cout << endl;
			cout << "All_reducing......" << endl;

			for(unsigned int i  = 0; i < _all_element.size(); i++)
			{
				_communicator.sum<Real>(_all_element[i]->_local_RD);
			}

			cout << "All_reduce completed" << endl;

			cout << endl;
			cout << "Normalizing......" << endl;

			_communicator.barrier();

			for (unsigned int ii=0;ii<_all_element.size();ii++)
			{
				for (unsigned int i=0;i<_all_element.size();i++)
				{
					_all_element[ii]->_local_RD[i] /= _particle_count;
				}
			}

			_communicator.barrier();

			cout << "Normalization completed" << endl;
		}

		else
		{
			ompComputeRD(_n_threads, _have_medium);
		}
		
		Deal_with_BSMC(_BSMC);

		cout << endl;
		cout << "Writing......" << endl;

		ofstream rd_write;
		rd_write.open(_filename.c_str(),ios::trunc);
		if(!rd_write.is_open())
			mooseError("Error opening file '" + _filename + "' fail.");

		rd_write << _all_element.size() << " " << _all_element.size() << endl;

		for (unsigned int ii=0;ii<_all_element.size();ii++)
		{
			for (unsigned int i=0;i<_all_element.size();i++)
			{
				rd_write << _all_element[ii]->_local_RD[i] << " " ;
			}
			rd_write << endl;
		}

		rd_write.close();

		cout << "Write completed" << endl;
	}

	else
	{
		cout << endl;
		cout << "Reading......" << endl;

		_communicator.barrier();

		for (unsigned int ii=0;ii<_all_element.size();ii++)
		{
			_all_element[ii]->_local_RD.resize(_all_element.size(), 0);
		}

		ifstream rd_read(_filename.c_str());
		if(!rd_read.good())
			mooseError("Error opening file '" + _filename + "' from  data.");

		string line;
		getline(rd_read, line);
		istringstream iss(line);
		int f;
		vector<int> head;
		while(iss >> f)
			head.push_back(f);

		unsigned int _num_i_of_RDij = head[0];
		unsigned int _num_j_of_RDij = head[1];

		if ( (_num_i_of_RDij!=_all_element.size()) || (_num_j_of_RDij!=_all_element.size()) )
			mooseError("Error file '" + _filename + "' connot march.");

		for (unsigned int i = 0; i < _all_element.size(); i++)
		{
			getline(rd_read, line);
			int j=0;
			istringstream transline(line);
			Real d;
			while (transline >> d)
			{
				_all_element[i]->_local_RD[j] = d;
				j+=1;
			}
		}

		_communicator.barrier();

		for(unsigned int i  = 0; i < _all_element.size(); i++)
		{
			_communicator.broadcast<Real>(_all_element[i]->_local_RD);
		}
		_communicator.barrier();

		cout << "Read completed" << endl;
	}

	cout << endl;
	cout << "/////////////////////Result of RD//////////////////////" << endl;

	for (unsigned int ii = 0; ii < _all_element.size(); ii++)
	{
		for (unsigned int i = 0; i < _all_element.size(); i++)
		{
			cout << _all_element[ii]->_local_RD[i] << " ";
		}
		cout << endl;
	}

	/*for (unsigned int i=0;i<_all_element.size();i++)
	{
		Real is_union = 0;
		cout << endl << "当前单元中心点：" << _all_element[i]->_elem->centroid() << endl;;
		cout << endl << "当前单元法向量：" << _all_element[i]->getSideElementNormal() << endl;;
		for (unsigned int j=0;j<_all_element.size();j++)
		{
			is_union += _all_element[i]->local_RD[j];
			cout << _all_element[j]->_elem->centroid() << "  RD:  " << _all_element[i]->local_RD[j] << endl;
		}
		cout << "is_union = " << is_union << endl;
	}*/

	end_time = clock();
	cout << endl;
	cout << "time: " << end_time-start_time << endl;

	if (_have_medium)
	{
		for (unsigned int i = 0; i < mymesh._userDefinedElem.size(); i++)
		{
			for (unsigned int j = 0; j < mymesh._userDefinedElem[i]._userDefinedSideElem.size(); j++)
			{
				delete mymesh._userDefinedElem[i]._userDefinedSideElem[j]._elem;
			}
		}
	}
}

void MCRadiationHeatFluxUserObject::initialize()
{
}

void MCRadiationHeatFluxUserObject::execute()
{
	int findi = Find_i(_current_side_elem);
	Real temp_pow4_bar(0);
	Real temp_pow3_bar(0);
	Real temp_bar(0);

	for(unsigned int _qp = 0; _qp < _q_point.size(); ++_qp)
	{
		temp_pow4_bar += (_JxW[_qp]*pow(_temperature[_qp],4) );
		temp_pow3_bar += (_JxW[_qp]*pow(_temperature[_qp],3) );
		temp_bar += (_JxW[_qp] * _temperature[_qp]);
	}

	temp_pow4_bar /= _current_side_volume;
	temp_pow3_bar /= _current_side_volume;
	temp_bar /= _current_side_volume;
	_temperature_pow4_bar[findi] = temp_pow4_bar;
	_temperature_pow3_bar[findi] = temp_pow3_bar;
	_temperature_bar[findi] = temp_bar;
}

void MCRadiationHeatFluxUserObject::finalize()
{
	{
		_communicator.sum<Real>(_temperature_pow4_bar);
		_communicator.sum<Real>(_temperature_pow3_bar);
		_communicator.sum<Real>(_temperature_bar);
	}

	for(int i = 0; i < _all_element.size(); i++)
	{
		_flux_radiation[i] = 0;
		_flux_radiation_jacobi[i] = 0;
	}

	computeRadiationFlux();

	{
		_communicator.sum<Real>(_flux_radiation);
		_communicator.sum<Real>(_flux_radiation_jacobi);
	}

	for(int i = 0; i < _all_element.size(); i++)
	{
		_temperature_pow4_bar[i] = 0;
		_temperature_pow3_bar[i] = 0;
		_temperature_bar[i] = 0;
	}
}

void MCRadiationHeatFluxUserObject::computeRD(int local_particle_count, bool have_medium)
{
	MooseRandom::seed( (2^24)*processor_id() );
	for(unsigned int ii  = 0; ii < _all_element.size(); ii++)
	{
		SideElement * cse = _all_element[ii];

		for (int j = 0; j < local_particle_count; j++)
		{
			int j_of_RDij = -1;

			j_of_RDij = Find_j_of_RDij(cse, _all_element, have_medium);
			
			if (j_of_RDij == -1)
				continue;
			else
				_all_element[ii]->_local_RD[j_of_RDij] = _all_element[ii]->_local_RD[j_of_RDij] + 1.0;
		}
		cout << "Computing " << ii << " of " << _all_element.size() << endl;
	}
}

void MCRadiationHeatFluxUserObject::ompComputeRD(int n_threads, bool have_medium)
{
#pragma omp parallel for num_threads(n_threads) \
	schedule(static,1)
	for(unsigned int ii = 0; ii < _all_element.size(); ii++)
	{
		for (unsigned int i = 0; i < _all_element.size(); i++)
		{
			_all_element[ii]->_local_RD[i] = 0;
		}

		SideElement * cse = _all_element[ii];

		for (int j = 0; j < _particle_count; j++)
		{
			int j_of_RDij = -1;

			j_of_RDij = Find_j_of_RDij(cse, _all_element, have_medium);
			
			if (j_of_RDij == -1)
				continue;

			else
				_all_element[ii]->_local_RD[j_of_RDij] = _all_element[ii]->_local_RD[j_of_RDij] + 1.0;
		}
		cout << "Computing " << ii << " of [" << _all_element.size() << "], " << 100.0 * ii / _all_element.size() << "% completed" << endl;
	}

	cout << endl;
	cout << "Normalizing......" << endl;
	for(unsigned int ii  = 0; ii < _all_element.size(); ii++)
	{
		for (unsigned int i = 0;i < _all_element.size();i++)
		{
			_all_element[ii]->_local_RD[i] = _all_element[ii]->_local_RD[i] / _particle_count;
		}
	}
	cout << "Normalization completed" << endl;
}



void MCRadiationHeatFluxUserObject::Deal_with_BSMC(int BSMC)
{
	for(int k = 0; k < BSMC; ++k)
	{
		vector<Real> N_i;
		N_i.resize(_all_element.size(), 1.0);
		Real temp = 0;
		for (unsigned int i=0;i<_all_element.size();i++)
		{
			Real normalized_parameter = _all_element[i]->_elem->volume()*_all_element[i]->_absorptivity;
			N_i[i] = N_i[i]*normalized_parameter;
			for (unsigned int j=0;j<_all_element.size();j++)
			{
				_all_element[i]->_local_RD[j] = _all_element[i]->_local_RD[j]*normalized_parameter;
			}
		}
		for (unsigned int i=0;i<_all_element.size();i++)
		{
			for (unsigned int j=0;j<_all_element.size();j++)
			{
				N_i[i] += _all_element[j]->_local_RD[i];
			}
		}
		for (unsigned int i=0;i<_all_element.size();i++)
		{
			temp = 0;
			for (unsigned int j=i;j<_all_element.size();j++)
			{
				temp = _all_element[i]->_local_RD[j] + _all_element[j]->_local_RD[i];
				_all_element[i]->_local_RD[j] = temp/N_i[i];
				_all_element[j]->_local_RD[i] = temp/N_i[j];
			}
		}
	}
}

void MCRadiationHeatFluxUserObject::computeRadiationFlux()
{
	Real sigma = 5.67e-8;
	Real Flux_Rad = 0.0;

	int local_begin, local_end, local_n;
	int quotient = _all_element.size()/n_processors();
	int remainder = _all_element.size()%n_processors();
	if(processor_id() < remainder)
	{
		local_n = quotient + 1;
		local_begin = processor_id()*local_n;
		local_end = local_begin + local_n;
	}
	else
	{
		local_n = quotient;
		local_begin = processor_id()*local_n + remainder;
		local_end = local_begin + local_n;
	}

	for (int i = local_begin; i < local_end; i++)
	{
		Flux_Rad = 0.0;
		for (unsigned int j = 0;j < _all_element.size(); j++)
		{
			Flux_Rad += (_all_element[j]->_local_RD[i])*(_all_element[j]->_elem->volume())*_all_element[j]->_absorptivity*_temperature_pow4_bar[j];
		}
		_flux_radiation[i] = sigma*Flux_Rad/_all_element[i]->_elem->volume();
		_flux_radiation_jacobi[i] = 4*sigma*(_all_element[i]->_local_RD[i])*_all_element[i]->_absorptivity*_temperature_pow3_bar[i];
	}
}

bool MCRadiationHeatFluxUserObject::ElemHaveNeighborInBlock(const Elem * elem, set<SubdomainID> block_ids)
{
	if (elem)
	{
		set<SubdomainID>::iterator it;
		for(it = block_ids.begin(); it != block_ids.end(); it++)
		{
			if(elem->subdomain_id() == *it)
			{
				return true;
			}

			else
			{
				continue;
			}
		}
		return false;
	}

	else
	{
		return false;
	}
}

bool MCRadiationHeatFluxUserObject::ElemInBlock(SubdomainID curr_SubdomainID, set<SubdomainID> block_ids)
{
	for (set<SubdomainID>::iterator it = block_ids.begin(); it != block_ids.end(); ++it)
	{
		if ( curr_SubdomainID != *it )
			continue;
		return true;
	}
	return false;
}

int MCRadiationHeatFluxUserObject::Find_i(const Elem * elem) const
{
	int findi = -1;

	for (unsigned int i = 0; i < _all_element.size(); i++)
	{
		if( (_all_element[i]->_elem->centroid() - elem->centroid()).norm() < RaylineZero)
		{
			findi = i;
			break;
		}
	}

	if(findi != -1)
		return findi;
	mooseError("Cannot find i.");
}

int MCRadiationHeatFluxUserObject::Which_SideelementIntersectedByLine(RayLine& ray, SideElement * sideelement_i, vector<SideElement*> sideelement_vec, Point & point)
{
	int j_max = sideelement_vec.size();
	int j_wanted = -1;
	Point pp = ray._p0;
	point = ray._p1;

	for(int j = 0; j < j_max; j++)
	{
		if( ((sideelement_vec[j]->getSideElementNormal()) * (ray._normal) > RaylineZero)  )
			continue;

		else if(!(ray.sideIntersectedByLine(sideelement_vec[j]->_elem, pp)))
			continue;

		else if((pp - ray._p0).norm() > (point - ray._p0).norm())
			continue;

		else
		{
			j_wanted = j;
			point = pp;
		}
	}
	return j_wanted;
}

int MCRadiationHeatFluxUserObject::Find_j_of_RDij(SideElement * sideelement_i, vector<SideElement*> sideelement_vec, bool have_medium)
{
	int j = 0;
	int j_of_RDij = -1;
	int k = 0;
	bool charge = true;
	SideElement * current_elem = sideelement_i;
	RayLine rayline_in;
	RayLine rayline_out;
	RayLine* ray_in = &rayline_in;
	Point p(0,0,0);
	rayline_in = (*current_elem).sendRay();
	while (charge && (k < _max_reflect_count) )
	{		
		if(have_medium)
		{
			j = findFinalSideId(rayline_in, p, current_elem);
		}
		else
		{
			j = Which_SideelementIntersectedByLine( rayline_in, current_elem, sideelement_vec, p);
		}
		
		double rand_num = MooseRandom::rand();

		if(j == -1 || rand_num <= sideelement_vec[j]->_transmissivity)
			return -1;

		else if(rand_num <= (sideelement_vec[j]->_absorptivity + sideelement_vec[j]->_transmissivity))
		{
			charge = false;
			j_of_RDij = j;
			break;
		}

		else if(MooseRandom::rand() <= sideelement_vec[j]->_diffuse_reflectivity)
		{
			rayline_out = sideelement_vec[j]->diffuseReflectRay(p);
			current_elem = sideelement_vec[j];
			rayline_in = rayline_out;
			j_of_RDij = j;
			k++;
			continue;
		}

		else
		{
			rayline_out = sideelement_vec[j]->mirrorsReflectRay(ray_in,p);
			current_elem = sideelement_vec[j];
			rayline_in = rayline_out;
			j_of_RDij = j;
			k++;
			continue;
		}
	}

	if(!charge)
		return j_of_RDij;

	else if(k == _max_reflect_count)
		return j_of_RDij;

	else
		return -1;
}

int MCRadiationHeatFluxUserObject::sideIntersectedByLine(const Elem * elem, int not_side, RayLine & ray_line, Point & intersection_point)
{
	unsigned int n_sides = elem->n_sides();

	// A Point to pass to the intersection method
	//  Point intersection_point;
	// Whether or not they intersect
	bool intersect = false;

	unsigned int dim = elem->dim();

	for (unsigned int i=0; i < n_sides; i++)
	{
		if (static_cast<int>(i) == not_side) // Don't search the "not_side"
			continue;

		// Get a simplified side element
		UniquePtr<const Elem> side_elem = elem->side_ptr(i);

		if (dim == 3)
		{
			// Make a plane out of the first three nodes on the side
			Plane plane(side_elem->point(0), side_elem->point(1), side_elem->point(2));

			// See if they intersect
			intersect = ray_line.intersect(plane, intersection_point);
		}
		else if (dim == 2)
		{
			// Make a Line Segment out of the first two nodes on the side
			RayLine ray_line_side(side_elem->point(0), side_elem->point(1), 1);

			// See if they intersect
			intersect = ray_line.intersect(ray_line_side, intersection_point);
		}
		else // 1D
		{
			// See if the line segment contains the point
			intersect = ray_line.contains_point(side_elem->point(0));

			// If it does then save off that one point as the intersection point
			if (intersect)
				intersection_point = side_elem->point(0);
		}

		if (intersect)
			if (side_elem->contains_point(intersection_point, RaylineZero))
				return i;
	}

// Didn't find one
	return -1;
}

/**
 * Returns the side number for elem that neighbor is on
 *
 * Returns -1 if the neighbor can't be found to be a neighbor
 */
int MCRadiationHeatFluxUserObject::sideNeighborIsOn(const UserDefinedElem * elem, UserDefinedElem * neighbor)
{
  unsigned int n_sides = neighbor->_elem->n_sides();

  for (unsigned int i = 0; i < n_sides; i++)
  {
	  if (neighbor->_userDefinedSideElem[i]._right_element == elem)
		  return i;
  }

  return -1;
}

int MCRadiationHeatFluxUserObject::pointInWhichSide(const Elem * elem, Point & point)
{
	unsigned int n_sides = elem->n_sides();

	for (unsigned int i = 0; i < n_sides; i++)
	{
		if( (*(elem->side_ptr(i))).contains_point(point) )
			return i;
	}
	return -1;
}

int MCRadiationHeatFluxUserObject::findFinalSideId(RayLine & ray_line, Point & point, SideElement * sideelement)
{
	int incoming_side = -1;
	UserDefinedElem * neighbor;

	const UserDefinedElem * first_elem = sideelement->_belong_to_which_elem;

	if (!first_elem)
	{
		cout << "first_elem" << endl;
		return -1;
	}

	else
	{
		const UserDefinedElem * current_elem = first_elem;
		incoming_side = sideelement->_is_which_sideelem;
		bool charge = true;
		while (charge)
		{
			int intersected_side = sideIntersectedByLine(current_elem->_elem, incoming_side, ray_line, point);

			if (intersected_side != -1) // -1 means that we didn't find any side
			{
				neighbor = current_elem->_userDefinedSideElem[intersected_side]._right_element;

				if (neighbor)
				{
					incoming_side = sideNeighborIsOn(current_elem, neighbor);
					current_elem = neighbor;
				}

				else
				{
					return current_elem->_haveWhichSideElement[intersected_side];
				}
			}
			else
			{
				/*如需增加介质，在此返回吸收单元*/
				return -1;
			}
		}
	}
	return -1;
}

void MCRadiationHeatFluxUserObject::print_mesh_to_file()
{
	string _filename = "mesh_for_fast_flux.dat";

	ofstream rd_write;
	rd_write.open(_filename.c_str(), ios::trunc);
	if(!rd_write.is_open())
    mooseError("Error opening file '" + _filename + "' fail.");

	rd_write << _all_element[0]->_elem->dim() << '	' << _all_element.size() << endl;
	for(unsigned int ii = 0; ii < _all_element.size(); ii++)
	{
		Point p = Point(0.0, 0.0, 0.0);
		int n_nodes = 0;
		if(_all_element[ii]->_elem->type() == TRI3)
		{
			n_nodes = 3;
		}
		else if(_all_element[ii]->_elem->type() == QUAD4)
		{
			n_nodes = 4;
		}
		else
		{
			mooseError("Error! This sideelement is not TRI3 or QUAD4. Please check it!");
		}
		rd_write << n_nodes << endl;
		for(unsigned int jj = 0; jj < n_nodes; jj++)
		{
			p = _all_element[ii]->_elem->point(jj);
			rd_write << p(0) << '	' << p(1) << '	' << p(2) << endl;
		}
		p = _all_element[ii]->getSideElementNormal();
		rd_write << p(0) << '	' << p(1) << '	' << p(2) << endl;
		rd_write << _all_element[ii]->_transmissivity << '	' << _all_element[ii]->_absorptivity << '	'
             << _all_element[ii]->_diffuse_reflectivity << '	'
             << _all_element[ii]->_mirrors_reflectivity << endl;
	}

	rd_write.close();

	cout << "Write elem_T_list to " << _filename << " completed" << endl;
}

void MCRadiationHeatFluxUserObject::print_T_to_file()
{
	string _filename = "T_for_fast_flux.dat";
	Point p = Point(0.0, 0.0, 0.0);

	ofstream rd_write;
	rd_write.open(_filename.c_str(), ios::trunc);
	if(!rd_write.is_open())
		mooseError("Error opening file '" + _filename + "' fail.");

	rd_write << _all_element.size() << endl;
	for(unsigned int i = 0; i < _all_element.size(); i++)
	{
		p = _all_element[i]->_elem->centroid();
		rd_write << p(0) << '	' << p(1) << '	' << p(2) << '	' << _temperature_bar[i] << endl;
	}

	rd_write.close();
}

void MCRadiationHeatFluxUserObject::print_flux_to_file()
{
	string _filename = "receive_flux.dat";
	Point p = Point(0.0, 0.0, 0.0);

	ofstream rd_write;
	rd_write.open(_filename.c_str(), ios::trunc);
	if(!rd_write.is_open())
		mooseError("Error opening file '" + _filename + "' fail.");

	rd_write << _all_element.size() << endl;
	for(unsigned int ii = 0; ii < _all_element.size(); ii++)
	{
		p = _all_element[ii]->_elem->centroid();
		rd_write << p(0) << '	' << p(1) << '	' << p(2) << '	' << _flux_radiation[ii] << endl;
	}

	rd_write.close();
}
